{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "identical-triangle",
   "metadata": {},
   "source": [
    "# High Molecular Weight Petroleum Pseudocomponents"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "responsible-pursuit",
   "metadata": {},
   "source": [
    "Thermo is a general phase equilibrium engine, and if the user provides enough properties for the components, there is no issue adding your own components. In this basic example below, a made-up extended gas analysis is used to specify a gas consisting of the standard real components and three heavier fractions, C10+, C12+ and C15+.\n",
    "\n",
    "A bare minimum of basic properties are estimated using the Kesler-Lee method (1976), and the estimated fraction molecular weights are turned into atomic compositions. The heat capacities of each pseudocomponent is found with\n",
    "the similarity variable concept of Lastovka and Shaw (2013) based on atomic composition.\n",
    "\n",
    "This example ends with calculating a flash at 270 Kelvin and 1 bar."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "excess-truck",
   "metadata": {},
   "outputs": [],
   "source": [
    "from math import log, exp\n",
    "import numpy as np\n",
    "from scipy.constants import psi\n",
    "from thermo import *\n",
    "from chemicals import *\n",
    "\n",
    "def Tc_Kesler_Lee_SG_Tb(SG, Tb):\n",
    "    r'''Estimates critical temperature of a hydrocarbon compound or petroleum\n",
    "    fraction using only its specific gravity and boiling point, from\n",
    "    [1]_ as presented in [2]_.\n",
    "\n",
    "    .. math::\n",
    "        T_c = 341.7 + 811.1SG + [0.4244 + 0.1174SG]T_b\n",
    "        + \\frac{[0.4669 - 3.26238SG]10^5}{T_b}\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    SG : float\n",
    "        Specific gravity of the fluid at 60 degrees Farenheight [-]\n",
    "    Tb : float\n",
    "        Boiling point the fluid [K]\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    Tc : float\n",
    "        Estimated critical temperature [K]\n",
    "\n",
    "    Notes\n",
    "    -----\n",
    "    Model shows predictions for Tc, Pc, MW, and omega.\n",
    "    Original units in degrees Rankine.\n",
    "\n",
    "    Examples\n",
    "    --------\n",
    "    Example 2.2 from [2]_, but with K instead of R.\n",
    "\n",
    "    >>> Tc_Kesler_Lee_SG_Tb(0.7365, 365.555)\n",
    "    545.0124354151242\n",
    "\n",
    "    References\n",
    "    ----------\n",
    "    .. [1] Kesler, M. G., and B. I. Lee. \"Improve Prediction of Enthalpy of\n",
    "       Fractions.\" Hydrocarbon Processing (March 1976): 153-158.\n",
    "    .. [2] Ahmed, Tarek H. Equations of State and PVT Analysis: Applications\n",
    "       for Improved Reservoir Modeling. Gulf Pub., 2007.\n",
    "    '''\n",
    "    Tb = 9/5.*Tb # K to R\n",
    "    Tc = 341.7 + 811.1*SG + (0.4244 + 0.1174*SG)*Tb + ((0.4669 - 3.26238*SG)*1E5)/Tb\n",
    "    Tc = 5/9.*Tc # R to K\n",
    "    return Tc\n",
    "\n",
    "def Pc_Kesler_Lee_SG_Tb(SG, Tb):\n",
    "    r'''Estimates critical pressure of a hydrocarbon compound or petroleum\n",
    "    fraction using only its specific gravity and boiling point, from\n",
    "    [1]_ as presented in [2]_.\n",
    "\n",
    "    .. math::\n",
    "        \\ln(P_c) = 8.3634 - \\frac{0.0566}{SG} - \\left[0.24244 + \\frac{2.2898}\n",
    "        {SG} + \\frac{0.11857}{SG^2}\\right]10^{-3}T_b\n",
    "        + \\left[1.4685 + \\frac{3.648}{SG} + \\frac{0.47227}{SG^2}\\right]\n",
    "        10^{-7}T_b^2-\\left[0.42019 + \\frac{1.6977}{SG^2}\\right]10^{-10}T_b^3\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    SG : float\n",
    "        Specific gravity of the fluid at 60 degrees Farenheight [-]\n",
    "    Tb : float\n",
    "        Boiling point the fluid [K]\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    Pc : float\n",
    "        Estimated critical pressure [Pa]\n",
    "\n",
    "    Notes\n",
    "    -----\n",
    "    Model shows predictions for Tc, Pc, MW, and omega.\n",
    "    Original units in degrees Rankine and psi.\n",
    "\n",
    "    Examples\n",
    "    --------\n",
    "    Example 2.2 from [2]_, but with K instead of R and Pa instead of psi.\n",
    "\n",
    "    >>> Pc_Kesler_Lee_SG_Tb(0.7365, 365.555)\n",
    "    3238323.346840464\n",
    "\n",
    "    References\n",
    "    ----------\n",
    "    .. [1] Kesler, M. G., and B. I. Lee. \"Improve Prediction of Enthalpy of\n",
    "       Fractions.\" Hydrocarbon Processing (March 1976): 153-158.\n",
    "    .. [2] Ahmed, Tarek H. Equations of State and PVT Analysis: Applications\n",
    "       for Improved Reservoir Modeling. Gulf Pub., 2007.\n",
    "    '''\n",
    "    Tb = 9/5.*Tb # K to R\n",
    "    Pc = exp(8.3634 - 0.0566/SG - (0.24244 + 2.2898/SG + 0.11857/SG**2)*1E-3*Tb\n",
    "    + (1.4685 + 3.648/SG + 0.47227/SG**2)*1E-7*Tb**2\n",
    "    -(0.42019 + 1.6977/SG**2)*1E-10*Tb**3)\n",
    "    Pc = Pc*psi\n",
    "    return Pc\n",
    "\n",
    "def MW_Kesler_Lee_SG_Tb(SG, Tb):\n",
    "    r'''Estimates molecular weight of a hydrocarbon compound or petroleum\n",
    "    fraction using only its specific gravity and boiling point, from\n",
    "    [1]_ as presented in [2]_.\n",
    "\n",
    "    .. math::\n",
    "        MW = -12272.6 + 9486.4SG + [4.6523 - 3.3287SG]T_b + [1-0.77084SG\n",
    "        - 0.02058SG^2]\\left[1.3437 - \\frac{720.79}{T_b}\\right]\\frac{10^7}{T_b}\n",
    "        + [1-0.80882SG + 0.02226SG^2][1.8828 - \\frac{181.98}{T_b}]\n",
    "        \\frac{10^{12}}{T_b^3}\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    SG : float\n",
    "        Specific gravity of the fluid at 60 degrees Farenheight [-]\n",
    "    Tb : float\n",
    "        Boiling point the fluid [K]\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    MW : float\n",
    "        Estimated molecular weight [g/mol]\n",
    "\n",
    "    Notes\n",
    "    -----\n",
    "    Model shows predictions for Tc, Pc, MW, and omega.\n",
    "    Original units in degrees Rankine.\n",
    "\n",
    "    Examples\n",
    "    --------\n",
    "    Example 2.2 from [2]_, but with K instead of R and Pa instead of psi.\n",
    "\n",
    "    >>> MW_Kesler_Lee_SG_Tb(0.7365, 365.555)\n",
    "    98.70887589833501\n",
    "\n",
    "    References\n",
    "    ----------\n",
    "    .. [1] Kesler, M. G., and B. I. Lee. \"Improve Prediction of Enthalpy of\n",
    "       Fractions.\" Hydrocarbon Processing (March 1976): 153-158.\n",
    "    .. [2] Ahmed, Tarek H. Equations of State and PVT Analysis: Applications\n",
    "       for Improved Reservoir Modeling. Gulf Pub., 2007.\n",
    "    '''\n",
    "    Tb = 9/5.*Tb # K to R\n",
    "    MW = (-12272.6 + 9486.4*SG + (4.6523 - 3.3287*SG)*Tb + (1.-0.77084*SG - 0.02058*SG**2)*\n",
    "    (1.3437 - 720.79/Tb)*1E7/Tb + (1.-0.80882*SG + 0.02226*SG**2)*\n",
    "    (1.8828 - 181.98/Tb)*1E12/Tb**3)\n",
    "    return MW\n",
    "\n",
    "def omega_Kesler_Lee_SG_Tb_Tc_Pc(SG, Tb, Tc=None, Pc=None):\n",
    "    r'''Estimates accentric factor of a hydrocarbon compound or petroleum\n",
    "    fraction using only its specific gravity and boiling point, from\n",
    "    [1]_ as presented in [2]_. If Tc and Pc are provided, the Kesler-Lee\n",
    "    routines for estimating them are not used.\n",
    "\n",
    "    For Tbr > 0.8:\n",
    "    .. math::\n",
    "        \\omega = -7.904 + 0.1352K - 0.007465K^2 + 8.359T_{br}\n",
    "        + ([1.408-0.01063K]/T_{br})\n",
    "\n",
    "    Otherwise:\n",
    "    .. math::\n",
    "        \\omega = \\frac{-\\ln\\frac{P_c}{14.7} - 5.92714 + \\frac{6.09648}{T_{br}}\n",
    "        + 1.28862\\ln T_{br} - 0.169347T_{br}^6}{15.2518 - \\frac{15.6875}{T_{br}}\n",
    "         - 13.4721\\ln T_{br} + 0.43577T_{br}^6}\n",
    "\n",
    "        K = \\frac{T_b^{1/3}}{SG}\n",
    "\n",
    "        T_{br} = \\frac{T_b}{T_c}\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    SG : float\n",
    "        Specific gravity of the fluid at 60 degrees Farenheight [-]\n",
    "    Tb : float\n",
    "        Boiling point the fluid [K]\n",
    "    Tc : float, optional\n",
    "        Estimated critical temperature [K]\n",
    "    Pc : float, optional\n",
    "        Estimated critical pressure [Pa]\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    omega : float\n",
    "        Acentric factor [-]\n",
    "\n",
    "    Notes\n",
    "    -----\n",
    "    Model shows predictions for Tc, Pc, MW, and omega.\n",
    "    Original units in degrees Rankine and psi.\n",
    "\n",
    "    Examples\n",
    "    --------\n",
    "    Example 2.2 from [2]_, but with K instead of R and Pa instead of psi.\n",
    "\n",
    "    >>> omega_Kesler_Lee_SG_Tb_Tc_Pc(0.7365, 365.555, 545.012, 3238323.)\n",
    "    0.306392118159797\n",
    "\n",
    "    References\n",
    "    ----------\n",
    "    .. [1] Kesler, M. G., and B. I. Lee. \"Improve Prediction of Enthalpy of\n",
    "       Fractions.\" Hydrocarbon Processing (March 1976): 153-158.\n",
    "    .. [2] Ahmed, Tarek H. Equations of State and PVT Analysis: Applications\n",
    "       for Improved Reservoir Modeling. Gulf Pub., 2007.\n",
    "    '''\n",
    "    if Tc is None:\n",
    "        Tc = Tc_Kesler_Lee_SG_Tb(SG, Tb)\n",
    "    if Pc is None:\n",
    "        Pc = Pc_Kesler_Lee_SG_Tb(SG, Tb)\n",
    "    Tb = 9/5.*Tb # K to R\n",
    "    Tc = 9/5.*Tc # K to R\n",
    "    K = Tb**(1/3.)/SG\n",
    "    Tbr = Tb/Tc\n",
    "    if Tbr > 0.8:\n",
    "        omega = -7.904 + 0.1352*K - 0.007465*K**2 + 8.359*Tbr + ((1.408-0.01063*K)/Tbr)\n",
    "    else:\n",
    "        omega = ((-log(Pc/101325.) - 5.92714 + 6.09648/Tbr + 1.28862*log(Tbr)\n",
    "        - 0.169347*Tbr**6) / (15.2518 - 15.6875/Tbr - 13.4721*log(Tbr) +0.43577*Tbr**6))\n",
    "    return omega\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "apart-cycling",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Basic composition and names. All pure component properties are obtained from Chemicals and Thermo.\n",
    "pure_constants = ChemicalConstantsPackage.constants_from_IDs(\n",
    "    ['water', 'hydrogen', 'helium', 'nitrogen', 'carbon dioxide', 'hydrogen sulfide', 'methane',\n",
    "    'ethane', 'propane', 'isobutane', 'n-butane', 'isopentane', 'n-pentane', 'hexane', \n",
    "    'heptane', 'octane', 'nonane'])\n",
    "pure_fractions = [.02, .00005, .00018, .009, .02, .002, .82, .08, .031,\n",
    "                  .009, .0035, .0033, .0003, .0007, .0004, .00005, .00002]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "another-footwear",
   "metadata": {},
   "outputs": [],
   "source": [
    "pseudo_names = ['C10-C11', 'C12-C14', 'C15+']\n",
    "pseudo_carbon_numbers = [10.35, 12.5, 16.9]\n",
    "pseudo_SGs = [.73, .76, .775] # Specific gravity values are based of the alkane series\n",
    "pseudo_Tbs = [447, 526, 589] \n",
    "\n",
    "# Using the estimation methods defined earlier, we obtain some critical properties\n",
    "pseudo_Tcs = [Tc_Kesler_Lee_SG_Tb(SG, Tb) for SG, Tb in zip(pseudo_SGs, pseudo_Tbs)]\n",
    "pseudo_Pcs = [Pc_Kesler_Lee_SG_Tb(SG, Tb) for SG, Tb in zip(pseudo_SGs, pseudo_Tbs)]\n",
    "pseudo_MWs = [MW_Kesler_Lee_SG_Tb(SG, Tb) for SG, Tb in zip(pseudo_SGs, pseudo_Tbs)]\n",
    "pseudo_omegas = [omega_Kesler_Lee_SG_Tb_Tc_Pc(SG, Tb) for SG, Tb in zip(pseudo_SGs, pseudo_Tbs)]\n",
    "# Estimate the hydroen counts\n",
    "hydrogen_counts = [(MW - C*periodic_table.C.MW)/periodic_table.H.MW \n",
    "                  for C, MW in zip(pseudo_carbon_numbers, pseudo_MWs)]\n",
    "# Get the atomic compositions\n",
    "pseudo_atoms = [{'C': C, 'H': H} for C, H in zip(pseudo_carbon_numbers, hydrogen_counts)]\n",
    "# Calculate the similarity variable of each species\n",
    "similarity_variables = [similarity_variable(atoms=atoms) for atoms in pseudo_atoms]\n",
    "\n",
    "pseudo_fractions = [.0003, .00015, .00005]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "restricted-puppy",
   "metadata": {},
   "outputs": [],
   "source": [
    "pseudos = ChemicalConstantsPackage(names=pseudo_names, MWs=pseudo_MWs, Tbs=pseudo_Tbs,\n",
    "                                   atomss=pseudo_atoms,\n",
    "                                  Tcs=pseudo_Tcs, Pcs=pseudo_Pcs, omegas=pseudo_omegas,\n",
    "                                  similarity_variables=similarity_variables)\n",
    "# Add the pure components and the pseudocomponents to create a new package of constant values\n",
    "# which will be used by the phase and flash objects\n",
    "constants = pure_constants + pseudos\n",
    "# Obtain the temperature and pressure dependent objects\n",
    "properties = PropertyCorrelationsPackage(constants=constants)\n",
    "# This is the feed composition\n",
    "zs = normalize(pure_fractions + pseudo_fractions)\n",
    "T = 270 # K\n",
    "P = 1e5 # bar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "continuing-plain",
   "metadata": {},
   "outputs": [],
   "source": [
    "kijs = np.zeros((constants.N, constants.N)).tolist() # kijs left as zero in this example\n",
    "eos_kwargs = dict(Tcs=constants.Tcs, Pcs=constants.Pcs, omegas=constants.omegas, kijs=kijs)\n",
    "\n",
    "# The API SRK equation of state is used, but other cubic equations of state can be uesd instead\n",
    "gas = CEOSGas(APISRKMIX, eos_kwargs, HeatCapacityGases=properties.HeatCapacityGases, T=T, P=P, zs=zs)\n",
    "liq = CEOSLiquid(APISRKMIX, eos_kwargs, HeatCapacityGases=properties.HeatCapacityGases, T=T, P=P, zs=zs)\n",
    "liq2 = CEOSLiquid(APISRKMIX, eos_kwargs, HeatCapacityGases=properties.HeatCapacityGases, T=T, P=P, zs=zs)\n",
    "phase_list = [gas, liq, liq]\n",
    "\n",
    "# Set up the three phase flash engine\n",
    "flashN = FlashVLN(constants, properties, liquids=[liq, liq2], gas=gas)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "worst-america",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3, 0.9827041561275567, [0.016838840039984482, 0.00045700383245886604])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Do the flash, and get some properties\n",
    "res = flashN.flash(T=T, P=P, zs=zs)\n",
    "res.phase_count, res.gas_beta, res.liquids_betas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "asian-scheme",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-1961.5089633224945,\n",
       " 1989.3915447041693,\n",
       " 19.675910651652533,\n",
       " 1.0015480015807005e-05,\n",
       " 0.027074133060773903)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.H(), res.Cp_mass(), res.MW(), res.gas.mu(), res.gas.k()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "legitimate-norfolk",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(769.2525386053419, 599.2086838769081)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.heaviest_liquid.rho_mass(), res.lightest_liquid.rho_mass()"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}